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Abstract: We formulate a model of relativistic fermions moving in two Euclidean dimen- 
sions based on a tight-binding model of graphene. The eigenvalue spectrum of the resulting 
Dirac operator is solved numerically in smooth U(l) gauge field backgrounds carrying an 
integer- valued topological charge Q, and it is demonstrated that the resulting number of 
zero-eigenvalue modes is in accord with the Atiyah-Singer index theorem applied to two 
continuum flavors. A bilinear but gauge non-invariant chirality operator appropriate for 
distinguishing the topological zero modes is identified. When this operator is used to cal- 
culate Q, it is found that the maximum topological charge capable of being measured in 
this fashion scales with the perimeter of the lattice. Some concluding remarks compare 
these results to what is known for staggered lattice fermions. 
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1. Introduction 

Following pioneering work in the 1980's there has been a recent revival of interest in 
lattice fermion formulations with a minimal flavor content of chirally-symmetric fermions 
|2], ^, ||, |||. Stimulated by the excitations known to obey a quasi-relativistic Dirac equation 
in graphene (a monolayer of carbon atoms arranged in a honeycomb lattice which has 
recently been realised experimentally), Creutz Q devised a four-dimensional Euclidean 
lattice action describing two species of massless chirally-invariant fermion, each centred at 
a special location ±p M in momentum space. Creutz's action is formulated on a hypercubic 
lattice, with tunable parameters enabling the magnitude of p^ to be controlled. Borigi Q 
soon developed the idea, finding an action with two flavors located this time at the origin 
and at (^-, £-). In each case the flavor content is the minimum consistent with the 

celebrated no-go theorem governing lattice fermion actions which are local, unitary, and 
chirally symmetric H. This has led to hopes that these formulations could form the basis 
for an inexpensive alternative to overlap fermions in realistic lattice QCD simulations with 
two light quark flavors. 

The feature of these actions is that in the vicinity of the special "Dirac points" the 
eigenvalues E of a suitably-defined Hamiltonian operator can be written as E{p fjL — p„) = 
±-ftr|z(p)|, where z(p) can be written as z^e 1 ® in two Euclidean dimensions and as zq + iz.a 
in four dimensions [|J. In 2d z maps an S 1 surrounding the Dirac point to a complex 
phase: in Ad the analogous mapping is from S 3 to a quaternionic space. In either case the 
solution of the fermion Hamiltonian engenders a non-trivial wrapping, implying that the 
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surface must enclose a zero, so that the only consistent result on shrinking it to a point 
is E = 0. Hence the existence of the Dirac points, and the linear nature of the dispersion 
E(j3u —pu) i n their immediate vicinity, is topologically stable. 1 One might therefore hope 
that the desirable properties listed above persist even once interactions with gauge degrees 
of freedom are introduced. 

A difficulty has been pointed out in M; the actions of 0, 3, while chirally symmetric, 
break important hypercubic and discrete symmetries such as parity and time-reversal, 
meaning that they are not protected against acquiring non-covariant counterterms such 
as eg. ■07^ through quantum corrections. This introduces a severe fine-tuning problem 
to any practical simulation programme based on the original actions. Recently Creutz 
has proposed a refinement of the action of || in which it is claimed such effects can be 
mitigated to the point where perhaps they are manageable [JjJ. 

The purpose of the current paper is to investigate the interaction of these minimal chiral 
fermions with gauge fields, not via a perturbative approach to radiative corrections ||, 
but rather in a non-perturbative manner via their response to a globally-defined topological 
charge. It is well-known that in a gauge background with integer-valued topological charge 
Q, the spectrum of the Dirac operator D obeys the Atiyah-Singer index theorem 

Q = n+-n_, (1.1) 

where n+(—) denotes the number of positive (negative) chirality zero modes of D. We 
have been inspired by the classic paper || of Smit and Vink, who studied the extent to 
which ([O]) is obeyed for both staggered and Wilson lattice fermions. Accordingly, we will 
investigate the response of minimal chiral fermions in 2d by calculating the spectrum of D 
in a background U(l) gauge field corresponding to a quantised homogeneous flux, which 
can be shown to carry Q / 0. Because of the combination of apparent simplicity with 
high symmetry, we have chosen to track closely the original connection with graphene and 
hence formulate our fermions on a honeycomb lattice. In Ad the analogous "hyperdiamond" 
lattice corresponds to a particular choice of parameters in 0, and has been explored further 
in || . Despite the aesthetic appeal of the honeycomb, it presents technical challenges which 
we feel are worth reporting in some detail. We will show that the index theorem fll.lD is 
indeed satisfied, and find a definition of chirality ipaaip capable of distinguishing between 
zero and non-zero modes. The drawback is that it is not gauge-invariant, and may not 
even be universal. 

The special properties of fermions hopping on a honeycomb lattice have been known 
in the literature for many years, ever since Semenoff fllOfl noted the relativistic nature 
of the dispersion around the Dirac points and solved the resulting Dirac equation in a 
homogeneous magnetic field in the long-wavelength limit to expose topologically-stable zero 



modes, confirming a result found in continuum field theory by Jackiw [11]. More recently, 
the topological aspects of honeycomb fermions interacting with an external magnetic field 
have been used to account for the unconventional quantum Hall effect in graphene [O] . The 



1 Because of the periodicity of the Brillouin Zone, the Dirac points must appear in pairs so that the 
overall wrapping vanishes ^). 
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Figure 1: Honeycomb lattice 



index theorem has also been used to analyse the effects of point defects in graphene [13]. In 
these papers the authors either worked in the continuum limit of the lattice model or worked 
in momentum space. In this work, by contrast, we will examine the applicability of the 
continuum results to fermions defined on the finite lattices relevant for QCD simulations; 
our concerns will be the approach to and recovery of the continuum limit predictions, and 
the technical issues associated with defining a chirality operator referrred to in the previous 
paragraph. 

The remainder of the paper is organised as follows. In Sec. [2] we specify the honey- 
comb lattice, define the lattice Dirac operator, and show that in the long-wavelength limit 
an action describing two continuum Dirac flavors is recovered. Sec. ||] reviews the index 
theorem, and outlines how U(l) configurations with Q / may be constructed on the 
honeycomb, and Sec. || presents the Dirac spectrum calculated both for free fermions on 
the honeycomb, and for continuum fermions on backgrounds with Q 7^ 0. In Sec. || we then 
present numerical results for the spectrum calculated on lattices up to size 100 x 100 with 
Q 7^ 0. The definition used for the chirality of a mode differs from the naive expectation 
based on free fermions. Nonetheless, we will show that both the spectrum and the index 
calculated on the basis of this chirality match analytic expectations provided Q is not too 
large; interestingly, the maximum value of Q for which continuum results are reproduced 
turns out to scale with the perimeter of the lattice. Our concluding remarks in Sec. ^ will 
contrast what we have found with what is known for staggered lattice fermions. Some 
technical details concerning the definition of the Fourier transform on a finite honeycomb 
lattice are postponed to an Appendix. 



2. Lattice action 

It is convenient to begin our presentation using a Hamiltonian devised for physical (ie. 
2 + 1-dimensional) graphene |0] . It is assumed that on each site r of a honeycomb lattice 
there is a mobile electron which may hop to a neighbouring site under the constraints of 
the Pauli Exclusion Principle. Electron spin may be ignored for now; the tight-binding 
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Hamiltonian is then 



3 

H = -tJ2 S b1 (r)U(r, Bi)a(r + *) + a+(r + ^(r, s < )6(r), (2.1) 

r£B i=l 

where t is a hopping parameter, Si(i = 1, 2, 3) are the three vectors along the links as shown 
in Fig. U 

Bi = (l,0)/, s 2 = (-i^)/, s 3 = (-I,-^)Z, (2.2) 

and / is the honeycomb bond length. The sites labelled A and B belong to inequivalent 
sublattices, on which the operators a) (a) and (b) respectively create (destroy) fermions. 
The variable U (r, Sj) is a U(l)-valued gauge connection emerging from the B site at r along 
Sj. The Hamiltonian ( |2.1[ ) is thus invariant under U(l) gauge transformations. 

To expose the relativistic nature of the low-energy excitation spectrum, define H$ = 
H[U = 1] and transform to momentum space (we will refine our definition of the Fourier 
transformation during the course of what follows) : 

H ° = Yl (Hk)a j {k)b{k) + $*(k)tf(k)a(k)^} (2.3) 

k 

with 



i = JL 



(2.4) 



Consider a basis of Fock states \k±) = (v^) -1 ^^) ± |0> where a|0) = 6|0) = 0, 
{a'(k), a(k')} = 5 2 (k — k'), {a, a} = {a^,a^} = etc. It is straightforward to see that 
{k±\Ho\k±) = ±($(A;) + <!>*( k)) = ±E(k), and hence that the spectrum is symmetric about 
zero. At half-filling (one electron per site) the Fermi energy is thus at E = 0. 

The dispersion relation E(k) vanishes not at k = 0, but at the six corners of the first 
Brillouin zone, which is also a hexagon but rotated by 90° with respect to the cells of 
Fig. |]. The corners of this hexagon lie at the centres of inequivalent equilateral triangles 
formed from reciprocal lattice points; the Brillouin zone thus contains two inequivalent 
Dirac points around which an effective low-energy description can be built, which we will 
take as K± = (0, ±^7^). Around these points we can expand: 

HK ± +p) = ±v F [ Py =F i Px ] + 0(p 2 ) (2.5) 

where the Fermi velocity vf = It is now possible, by defining field operators in the 
neighbourhood of the Dirac points via a±(p) = a(K± + p) etc, to recast the Hamiltonian 
in relativistic form: 



v F ^2^(p)a.p^(p) (2.6) 
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where \I/ is the column vector a+, a_ , b-) T and the 4x4 matices a are defined by 

* = n )> -J 

so that {aj,<x,-} = 2(5^-. In this form i?o is easily seen to be proportional to the Dirac 
Hamiltonian describing a single massless four-component spinor moving with speed vp. 
For physical graphene, the Hamiltonian fl2,l|) must be modified to incorporate electron 
spin; this results in a relativistic d = 2 + 1 model with two four-component flavors. 

Starting from the Hamiltonian (|2.1| ) an action for chiral gauge theory in d = 2 + 1 was 



proposed by Jackiw and Pi in [ 14 \ . In this paper we instead recast it as a d = 2 Euclidean 
quantum field theory with action of the form S = X-Dx, describing two species of fermion 
field sitting on a honeycomb lattice, each species occupying a distinct sublattice. The 
resulting equation of motion resembles the Dirac equation in the long wavelength limit; 
the Dirac operator D can be written as 

(D X )(x) = Di{x + 6) X (x + 6) + D 1 (x- 6) X (x - 6) + D 2 (x + 2) X (x + 2) 

+D 2 (x-2)x(x-2) + D 3 (x) X (x). (2.8) 

In writing the operator this way we have introduced the notion of a lattice of identical A 
sites, with rhombus-shaped unit cells of side a = V3l each containing one A site and one B 
site (see Fig. ^). Each cell is indexed by a vector x, whose form will be specified below, and 
primitive vectors = si — S3 and 1 = S2 — si define the lattice axes. It is also convenient 
to define the dependent vector 2 = + 1. The different elements of the operator are then 
written 

( 0\ , fO U*(x-0,sx) 

V ; \U(x,s 2 ) OJ K V° 

D3{X) = U*,-s) J' (2 ' 



V(')- ( XA S ) (2-10) 



with the spinors 

where xa(x) and xb(^) are single-component Grassmann fields located at sites A and B 
of cell x respectively. 

Note that with our definition of the link directions all terms in the action are either of 
the form xbUxa or xaU*xb- While the antihermiticity of the Dirac operator ( [2.8| , |2.9| ) is 
not manifest, a little care and close inspection of Figs. Qj2| will convince the reader that this 
is indeed the case. The above analysis expanding S about the Dirac points k = K± goes 
through as before; writing ipi = (xb+,XA+) T , $1 = {xb+,XA+), ^2 = (~XA-, ~Xb-) T 
and $2 = {xa-,Xb-), with XA±(j>) = Xa{K± + p) as before, we obtain 

o, 2 

p ct=l 
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that is, a relativistically covariant action describing two flavors of two-component spinor 
moving in d = 2 Euclidean dimensions, each flavor localised at one single Dirac point. 
Chirality is then naively defined by the bilinear 

i>a<73i>a = XB+XB+ ~ XA+XA+ ~ XA-XA- + XB-XB- = XBXB ~ XAXA, (2.12) 

where the second equality assumes that all parts of momentum space can be treated uni- 
formly, leading in effect to a staggered order parameter. As we shall see in Sec. |5] below, 
the definition needs to be modified in the presence of gauge fields. 

To define a finite, translationally-invariant lattice we need to close the manifold by 
specifying boundary conditions. In this paper we have studied two distinct possibilities. 
The technically simpler choice is to close the manifold along the two non-orthogonal axes 
and 1, for instance defining an Lq x L\ system by requiring /(x + LqO) = f{x + L\l) = f(x). 
We call this the "primitive" lattice shown in red in Fig. ||. Note it contains LqLi distinct 
hexagons. We can alternatively choose to implement the boundary conditions along the 
orthogonal axes X and Y, in this case calling it the "perpendicular" lattice shown in black 
in Fig. [3|. Some care is needed in indexing the lattice this way: it is convenient to assign 
the two A sites r and r — S3 + S2 and the B site r — S3 the same X index, but to assign 
them Y indices of respectively eg. 0, ^, and 1. In this way a lattice which extends Lx 
units along X and Ly units along Y contains LxLy distinct hexagons. 

3. Index theorem 

Let us consider fermions in the presence of a topological charge Q created by a background 
gauge field configuration. In the continuum, the Atiyah-Singer index theorem relates the 
topological charge to the number of chiral zero modes of the fermion. If tpi are the eigen- 
states of the antihermitian Dirac operator with eigenvalue iEi, then 

Q= ^2 Vila^i = n+ -n-, (3.1) 

where is the number of zero eigenvalue modes of positive (negative) chirality, ie. 

satisfying 75^ = iV'i- A heuristic derivation of this relation for lattice fermions is given 




A 

Figure 2: The primitive cell 
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Figure 3: Perpendicular compactification 



Figure 4: Primitive compactification 



in 0. Since chiral symmetry is minimally broken, we anticipate the above relation holds 
good on a honeycomb lattice with a suitable generalisation of the chirality operator 75. 

In two Euclidean dimensions topological charge density is proportional to the magnetic 
field strength tensor F12 . The two-dimensional analogues of instantons are localised vortices 
carrying a quantised magnetic flux; the topological charge Q is defined by 



Q 



1 

2^ 



d 2 xF 12 . 



(3.2) 



On a finite system it is also possible to define homogeneous backgrounds with Q ^ 0. 
Consider the abelian field strength tensor F12 = d x A y — d y A x on a 2d system with bound- 
aries closed in orthogonal directions. Following ||, we then choose A x (x,y) = —coy and 
A y (x, y) = 0, so that F12 = uj. On a L x x L y square lattice of spacing a the gauge field A\ 
at the boundary y = L y a is related to that at y = by a gauge transformation S 



A x (y = 0) = A x (y 



L y a) + iVL y d x VL y 



(3.3) 



where 



n y (x,y) = e i ^ ax . 



(3.4) 



The discontinuity is permitted since any gauge invariant object remains continuous across 
the boundary of the lattice. Demanding periodicity in the x-direction as well results in a 
field strength quantised consistently with ( |3.2[ ): 



2tt 



L x L y a? 



Q. 



(3.5) 



On the honeycomb lattice the gauge background with constant field strength u over the 
lattice, and accordingly equal flux ojA through each hexagonal plaquette of area A 
is quantised according to 



47T 



uJ 



\/3L x L y a 2 



Q. 



(3.6) 
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The link field configuration depends on which boundary condition we consider. For the 
primitive lattice boundary condition, a possible choice is: 

U(x,si) = exp ( -i^-ujxoa 2 ] ; j7(a;,s 2 ) = 1; U(x, s 3 ) = 1, (3.7) 



2 

for all cells except those of the last row with xq = Lq — 1 where in addition we require 

U (x = L - 1, xx, s 3 ) = exp ^-i^o;LoXia 2 ^ . (3.8) 

It is readily checked that each hexagonal plaquette then has the value exp(iujA). 

For the perpendicular boundary condition, a link field configuration for the same con- 
stant field strength uj could be chosen as follows. For xy = 0, |, 1, . . . , Ly — |: 

U(x,s\) = exp ^—i^-ujxya 2 ^j ; U(x, S2) = 1; ?7(x,S3) = 1, (3.9) 
whereas for the links in the last row we need 

U(x Y = Ly - i,xx,s 2 ) = exp ^i^-ujLyxxa 2 ^ ; 

U(xy = Ly,xx,sz) = exp ^i^-ujLyXxo 2 ^ ■ (3.10) 

4. Dirac spectrum: analytic results 

With the above gauge field configuration we need to solve the lattice Dirac equation 

D X i = iEiXi (4.1) 

with Xi( x ) satisfying primitive or perpendicular boundary conditions as appropriate. For 
an arbitrary gauge field background we can do only that numerically. In this section we 
first discuss the spectrum on a honeycomb lattice without any gauge field, and then for a 
continuum Dirac operator on a smooth background field with Q 7^ of the kind outlined 
in Sec. |j| 

If we consider the boundary condition on the primitive lattice and in the free field 
limit, a plane wave ansatz gives the eigenvalues 



E(k ,h) = ±^3 + 2 cos k + 2cos(A;o - h) + 2 cos k[. (4.2) 

For periodic boundary conditions along and 1, the allowed momentum modes are given 
by 

7 2vrn 

k = — — , n = 0, 1,2, • • • L - 1; 
^0 
27rni 

h = -pA ni = 0,l,2,...Li-l. (4.3) 
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Note that koO and k\l are not orthogonal. One can clearly see that E(ko,k{) is not 
minimised at k$ = k\ = 0, but rather at the Dirac points K± = (^^,±^): recovering 
the continuum relativistic dispersion relation is therefore nontrivial. Writing k = K± + q, 
the leading term in the expansion is 



E(q , qi ) = ±^Jql - q qi + q\ . (4.4) 

Replacing the non-orthogonal qo, q± by orthogonal momenta p x , p y via f/o,i =Pn± "^Py 
we get the desired relativistic dispersion relation 



E(p x ,P y ) = ±^JPx 2 +P y 2 - (4.5) 

If expressions (4.2-4~5|) are required in terms of dimensionful momenta k,q,p, then note 
that p = ap etc. 

The solution for a single flavor of Dirac fermion moving in a 2d uniform background 
magnetic flux density u was first discussed in [11]; here we review the explict solution given 
in p. The equation is 

written 

2 

D ^3 = D H a ^3 = iE 3^3i ( 4 - 6 ) 
ft=l 

where tp is a two-component spinor. For a background with topological charge Q ^ there 
are \Q\ independent solutions of the form (with a = 1): 

Mx,y) ocf; e^A^^e-iM^^y+flOIMW ^\(y±^( j+ i\Q\))) ± (4.7) 

£=-oo V M 7 

where j = 0,1, . . . ,\Q\ — 1, 0+ = ( ^ J ' = f o) ' ^™ 8X6 ^" erm ^ e P or y nom i a l s °f order n, 



and Q and u are related via ( |3.2| ). The corresponding eigenvalues are given by 

^ = (2n+l)HTw. (4.8) 
Rearranging, we find a spectrum 

£& = 2m|w|, m = 0,1,2,... (4.9) 

with degeneracy 

^ = \2|g| ™>o. (410) 

The \Q\ zero modes are all proportional to (/>+ (</>-) for Q positive (negative), in accordance 
with the index theorem (3.1). For m > an equal number of positive and negative chirality 
solutions can be found. The increase of g m with uj is a relativistic analogue of the Landau 
levels observed in metals in a strong magnetic field. For the two continuum flavors described 
by the honeycomb Dirac operator fl2.§|), the index theorem thus predicts 2\Q\ zero modes, 
a result first obtained in [10]. 
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Figure 5: Eigenvalue spectrum as a function of Q on a 30 x 30 lattice (the vertical scale is measured 
in units where a~ 2 = -^). 



5. Numerical Results 



In order to analyze the spectrum of the Dirac operator in various gauge field backgrounds, 
the matrix —D 2 [U] was diagonalised via a subspace iteration technique, using Chebyshev 
polynomial iteration to accelerate the convergence of the eigenvalues E 2 . Since small 
eigenvalues converge at a faster rate than the high lying eigenvalues, locking the already 
converged eigenvalues and eigenvectors also accelerates the convergence of the other eigen- 
values. The locked eigenspaces are only used to orthogonalise the remaining subspaces. 
This algorithm is also suitable to find the few lowest lying eigenvectors. Further details 



may be found in [15, 16 1. 

In Fig. [5] we plot the 60 smallest eigenvalues E 2 calculated on a 30 x 30 lattice with 
primitive boundary conditions, for values of topological charge Q varying between 1 and 
26. The spectrum for the perpendicular boundary conditions is identical. On this lattice 
Lua 2 = 0.00698Q. Close inspection of the figure reveals very good agreement with both the 
eigenvalue prediction (|4.9| ) and the degeneracy pattern ( f4.10| ), recalling that for two flavors 
we expect a degeneracy 2g m . In particular, the triangular "carpet" with E 2 = corresponds 
to the zero modes with degeneracy growing linearly with Q as predicted by the index 
theorem ( |3.1[ ). It is important to note that for these smooth background configurations 
the zero-mode eigenvalue is equal to zero within machine precision, just as is the case for 
staggered fermions |8| . In what follows we will strengthen this correspondence by specifying 
a chirality operator appropriate for honeycomb fermions. 

In order to proceed, recall the discussion of Section and in particular that the low 
energy modes are located in the neighbourhood of the Dirac points, ie. away from the 
origin of momentum space. It is therefore helpful to perform analysis in momentum space 
by Fourier transforming the eigenvectors Xi( x )- K turns out that calculating the discrete 
Fourier transform on a honeycomb lattice is rather tricky ]l7||, essentially because the range 
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(a) E — 0, Multi-valued FT, A sublattice 
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(c) E ^ 0, Multi-valued FT, A sublattice 



(d) E ^0, Single valued FT, A sublattice 



Figure 6: Eigenvector profiles \\\ 2 plotted in momentum space. Results taken on a 30 x 30 lattice 
with primitive boundary conditions. 

of k- values required for a unique invertible Fourier transform to exist is larger than is the 
case for a square lattice. As shown in the Appendix, we have the choice of defining a single- 
valued transform ranging over either QLxLy (perpendicular) or 9LqL\ (primitive) modes, 
or a multi- valued transform ranging over 2LxLy (perpendicular) or LqL\ (primitive). This 
is exemplified in Fig. [f| where the single- valued transform shown in Fig. |d to a very good 
approximation consists of nine copies of the multi- valued transform defined over a smaller 
range shown in Fig. ^c. In fact, while we have carried out all subsequent analysis using 
both variants of the Fourier transform, the results in all cases are found to be identical, as 
exemplified by Fig. || below. 

In Fig. | we have chosen to transform the eigenvector densities on A and B sublattices 
separately to expose an important distinction between zero and non-zero modes. In each 
case the modes are tightly localised around two complementary locations in /c-space, which 
we identify with the ± Dirac points discussed previously. However, while the non-zero mode 
of Figs. |6|c,d clearly has support at both Dirac points, for the zero mode the eigenvector on 
the A sublattice is supported only near the + point (Fig. |6|a) whereas on the B sublattice 
it is supported only near the - point (Fig. ^b). 

This behaviour can be understood via the following heuristic argument. Consider 



solving the continuum problem on an infinite volume in a more symmetric gauge: A x 
— jy, A y = ijX. Solutions can readily be found of the form 



ij)e± oc (x t iy) exp f — —{x +y)\<t>± (5.1) 

with corresponding eigenvalues 

4 : = (H1)[Htw] i £ = 0,1,2,... (5.2) 



The spectrum (4.9 , 4. ld| ) is reproduced, with both chiralities contributing to non-zero 



modes, but with only positive chirality zero-modes present for u > and vice-versa. How- 
ever, in this gauge ipg± is also an eigenstate of orbital angular momentum L = i(yd x — xd y ) 
with eigenvalue ^f£. Physically, the describe particles executing circular motion (the 
modes ( |5.1| ) are localised on annular regions centred at the origin) with opposite senses for 
+ and — states - in other words the particle's charge and hence its response to a magnetic 
field is determined by its chirality. Now, in the rest frame both chiralities yield orbits 
of the same shape. However, for our honeycomb fermions the states are located at the 
Dirac points, and hence the previous picture needs to be Lorentz-boosted. A charged par- 
ticle moving through a magnetic field with non-zero linear momentum has as trajectory 
a 2d projection of a helix; + and — particle states can no longer be superimposed and 
must therefore be described by different wavefunctions. Hence a state with a well-defined 
chirality is necessarily localised around a single point in fc-space. 

We learn from this argument that constructing states with well-defined continuum 
quantum numbers may not be straightforward for honeycomb fermions. In particular, the 
naive definition of chirality ( |2.12 ) derived for free fermions appears not to be suitable for 



the fermion modes in constant background flux of Fig. ||a,b , since in this case it receives 
cancelling contributions from A and B sublattices. Instead, we propose the following 
definition of chirality for interacting honeycomb fermions: 

tpa^^Ja = -XB+XB+ + XA+XA+ ~ XA-XA- + XB-XB-- (5.3) 



Equation (5^3) is evaluated in Fourier space, with each mode counted as + or — depending 
on which Dirac point it lies closer to. 2 We note in passing that the chirality operator 



introduced in the (2+l)-dimensional treatment of 14] also assigns opposite chiral charges 



to fields localised around the different Dirac points. In Fig. [7| we plot the expectation 
values of the chirality operator (|5.3|) evaluated on the lowest 40 eigenmodes of a 30x30 
system with Q = 4. We have repeated the calculation, each time implementing a random 
gauge transformation of the form e l9 ^ p at each site, where 9 is uniformly distributed 
around the circle and p is a parameter. This transformation of course leaves the spectrum 
unchanged. For the smooth untransformed background the chirality (S3) = +1 to a good 
approximation for the 2\Q\ topological zero modes. For non-zero modes the chirality has 



2 The location of the Dirac point depends on the boundary conditions and in practice is determined by 
the maximum of \\- 
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Figure 7: Chirality vs. mode number for various random gauge noise p 



a smaller magnitude and a fluctuating sign; moreover its sum over all degenerate non- 
zero modes is exactly zero. Both of these are of course minimum requirements for a 
realistic chirality operator. However, Fig. [7] also confirms that the operator is not 

gauge invariant, which is not surprising since it is formulated in momentum space. As the 
amplitude of the short- wavelength noise injected into the gauge background grows with p, 
the magnitude of (S3) falls steadily, until eventually the zero and non-zero modes become 
indistinguishable. Using the definition (5.3) on backgrounds with p = we present results 



for the topological charge Q index , as evaluated via the index theorem (3.1) on the zero 



modes, versus the charge Q flux , obtained by integrating the background flux ( |3.2| ), for a 
range of lattice sizes using both boundary conditions in Fig. |8| For Q not too large the 
curves fall on a straight line of unit slope independent of lattice volume, confirming the 
validity of the definition (|5.3|) and verifying the index theorem. However beyond some 
value of Q, which depends on L, the curves reach a maximum and then fall with increasing 
flux density. Different behaviour is observed for the two kinds of boundary condition but 
reassuringly, as mentioned above, the results are insensitive to which definition of the 
Fourier transform is used. 

In fact, the maximal topological charge achievable on a given lattice depends linearly on 
the length of its perimeter, rather than the area as naively one would expect. For instance 
in the case of an Lx x Ly perpendicular compactification, provided that that a different 
scale factor for the X and Y directions is chosen, then the maximum achievable Q mdex is 
linearly proportional to both Lx and Ly independently. In particular in Fig. [9| we plot 
Qmax 1 as a function of "perimeter" P = 3Lx + ^Ly for all possible combinations of Lx 
and Ly in the range 20-30 independently, plus for other larger lattices with Lx = Ly. A 
possible interpretation of this phenomenon is that for a 2d U(l) gauge theory with constant 
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Figure 8: Q lndcx versus Q flux for different lattices 




magnetic flux and fixed Q, it is always possible to perform a gauge transformation that 
moves all information about the topological charge in the system to the border of the 



lattice, where it will be encoded by transition functions Q x , £l y such as (^4). In this case 
the topological charge can be identified with the number of windings of a scalar field around 



the border, consistent with the quantisation condition (3.5). It follows that the maximal 
resolution obtainable is given by the number of points on the perimeter, i.e. the ratio of 
the length of the perimeter to the step size, since we are not able to probe a field winding a 
greater number of times than the number of points defining the discrete Fourier transform. 



6. Concluding Remarks 

In this paper we have demonstrated that the spectrum of a simple fermion model formu- 
lated on a 2d honeycomb lattice does indeed reproduce one of the most important non- 
perturbative features of relativistic chiral fermions interacting with a gauge field, namely 
the index theorem ( |1.1| ) relating the number of zero-eigenvalue chiral modes to the back- 
ground topological charge. We have done this first by calculating the gauge-invariant 
spectrum in a particularly smooth gauge background and showing that it coincides with 
analytic expectations yielding zero modes with the correct degeneracy, and next by iden- 
tifying a suitable chirality operator S3 enabling the index to be calculated. The operator 
fl5.3| ) distinguishes between fermion fields located on differing A and B sublattices in real 
space, and between fields localised at differing ± Dirac points in momentum space. As such 
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Figure 9: Q^tT vs - perimeter P in the perpendicular case 



it is necessarily not gauge- invariant, and hence has limited practical value except in the 
artificially-constructed smooth backgrounds used here. Since the spectrum and by exten- 
sion det-D are gauge-invariant, however, this need not deal a fatal blow to any simulation 
programme based on honeycomb fermions. 

Since we have constructed the operator (|5^) to work on a smooth gauge background, 
it is legitimate to ask how universal it is, ie. how would it respond in a non-uniform 
distribution of topological charge, such as that found in the vicinity of a vortex line? In 
the deep continuum limit all gauge backgrounds can be reduced at least locally to a smooth 
background by a suitable gauge transformation, where we know (|5.3| ) is appropriate, so it 
is difficult to see how any alternative definition could be preferred. Nonetheless, it seems 
likely that (£3} as defined by ( |5.3D will be extremely susceptible to lattice artifacts. 

We have also studied the maximum Q 1 ^^ observable in a lattice simulation, and shown 
that this scales with the lattice perimeter, arguing that the limit is related to the maximum 
resolution of the discrete Fourier transform along a lattice direction. To our knowledge this 
is a new observation; it would be interesting to repeat this analysis for staggered lattice 
fermions on a square lattice. 

We close by contrasting what we have found to what is known for staggered lat- 
tice fermions. Superficially the two approaches are very similar; one starts with a single- 
component Grassmann field on each site and then reallocates the degrees of freedom into a 
new basis to recover fields with spin and flavor quantum numbers appropriate to continuum 
fermions. The method used by Smit and Vink [|j partitions the original field in momentum 
space, according to a formalism originally developed in |p^| . In that case, however, the 2 d 
"Dirac points" include the origin, and in contrast to the discussion leading to eqn. ( 2,11| ), 
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each continuum flavor is evenly distributed across all such points. Bilinears with a specified 
spin/flavor structure constructed from free fermion fields continue to be applicable in the 
presence of gauge fields. Moreover a distinct repartition of the fields, based on their relative 
location in local "hypercubes" each containing 2 d lattice points, is also possible [19], en- 
abling the construction of local gauge-invariant bilinears by suitable insertion of products 
of link variables. For instance the chirality operator for 2d staggered fermions couples \ 
and x & t opposite corners of an elementary square || . The two formalisms coincide in the 
long-wavelength limit |2C]. For honeycomb fermions, the continuum flavors are localised at 
different Dirac points. We have argued in Sec. || that this implies that bilinears appropriate 
for free fermions need not continue to be correct once gauge interactions are introduced; it is 
also the case that construction of bilinears such as (|3l|) requires a simultaneous reshuffling 
both in position and in momentum space, making gauge invariance impossible to achieve. 
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A. Fourier transform on a honeycomb lattice 
A.l Single- valued Fourier transform 

The main result we want to obtain in this section is to write down a unique Fourier 
transform function f(u,v) with a well-defined inverse. Given the peculiar structure of the 
hexagonal lattice, it will be shown that the range of values of u and v of f(u,v) is larger 
than for the square case, and depends on the compactification. 

To define a Fourier transform we need to expose which elementary translations leave 
the lattice unchanged; naturally the different types of compactification implement these 
translations differently. If we want a natural coordinate to enumerate the sites of the lattice 
we need to define a different support function for each compactification type. 

We can then define a support function that defines the location of all the sites of our 
lattice Tl(x,y), by dividing the sites into classes, where by a class we mean the subgroup 
of sites related by a integer combination of elementary translations. It's straightforward to 
see that for the perpendicular lattice we need four classes to reproduce the entire lattice: 

n_i_(x,2/) = ^2 5(x - ^ - 3n, y - V3m) + S(x - 3n, y - ^ - VHm) + 

n=0,L x -lm=0,L y /2-l 

1 \/3 
6(x — — — 3n, y — Vom) + 6(x — 2 — 3n, y v3m), (A.l) 

while for the primitive lattice we need only two classes: 

U p (x,y) = 5(x - V3n, y - V3m) + 

n=0,L -l n=0,Li-l 

5(x - V3n - ^, y - V3n - -U (A.2) 
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The Fourier transform of any function f(x,y) can then be defined as 
f(u,v) = yj dxdy U(x,y) f(x,y) e^, 
where p is again different for the two compactifications 

Perpendicular — ► p = (— — u X, -7=7 — v Y) 



(A.3) 
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This leads to the two different definitions of Fourier transform: 
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V3 

Now note that the periodicity of the transformed functions so defined differ from the 
original: 

h (0,0) = L(6L X ,L Y ) 

/ P (0,0) =/ p (3Lo,3Li). (A.7) 

To define the inverse transform we cannot define an unique formula for all classes. Rather, 
the inverse must be calculated on each class separately. 

Perpendicular: 



/± (- + 3n,v / 3m) = 7 _= 
/_!_(- + 3n, — V3m) = ... 
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A. 2 Multi-valued Fourier transform 



If we relax the requirement of having a single- valued function of momentum, but still wish 
to expose the different structure of the A and B sites in the Fourier transform, we obtain 
two further definitions. In this case the the range of values of u and v depends on the 
compactification, but is always smaller than the previous case. The procedure to obtain 
these function is totally equivalent to the previous case, so here we only show the results. 
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Perpendicular Inverse: 
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Primitive Inverse: 

Lo-lLl-l -(nu,mv\ 

f p (V3n,V3m) = -±= £ E e ^ ^ 



-7T=r E E e V Ll hi B) {u,v)e \ 3L ° S W 

(A.13) 



- 19 - 



References 

[1] L.H. Karsten, Phys. Lett. B 104, 315 (1981); 

F. Wilczek, Phys. Rev. Lett. 59, 2397 (1987). 

M. Creutz, JHEP 0804, 017 (2008). 
A. Borici, Phys. Rev. D 78, 074504 (2008). 

P.F. Bcdaque, M.I. Buchoff, B.C. Tiburzi and A. Walker-Loud, Phys. Lett. B 662, 449 
(2008). 

K. Cichy, J. Gonzalez Lopez, K. Jansen, A. Kujawa and A. Shindler, Nucl. Phys. B 800 
(2008) 94. 

H.B. Nielsen and M. Ninomiya, Nucl. Phys. B 185, 20 (1981) [Erratum-ibid. B 195, 541 
(1982)]; Nucl. Phys. B 193, 173 (1981). 

M. Creutz, arXiv:0808.0014 [hep-lat]. 

J. Smit and J. C. Vink, Nucl. Phys. B 286, 485 (1987). 

P.F. Bcdaque, M.I. Buchoff, B.C. Tiburzi and A. Walker-Loud, Phys. Rev. D 78, 017502 
(2008). 

G. W. Semenoff, Phys. Rev. Lett. 53 (1984) 2449. 

R. Jackiw, Phys. Rev. D 29 (1984) 2375 [Erratum-ibid. D 33 (1986) 2500]. 

V.P. Gusynin and S.G. Sharapov, Phys. Rev. Lett. 95 (2005) 146801; 
Y. Hatsugai, T. Fukui and H. Suzuki, Phys. Rev. B 74, 205414 (2006). 

J.K. Pachos and M. Stone, Int. J. Mod. Phys. B21 (2007) 5113. 

R. Jackiw and S. Y. Pi, Phys. Rev. Lett. 98, 266402 (2008). 

Y. Saad, Numerical Methods For Large Eigenvalue Problems, Manchester University Press, 
1992. 

L. Del Dcbbio, L. Giusti, M. Liischcr, R. Pctronzio and N. Tantalo, JHEP 0602, 011 (2006). 

U.R. Alim and T. Moller, A discrete Fourier transform for the hexagonal and body-centered 
lattices, SFU Computing Science Technical Report 2008-14; 
A. Vincc and X. Zheng, J. Math. Imaging Vis. 28, 125 (2007). 

C. van den Doel and J. Smit, Nucl. Phys. B 228, 122 (1983). 

H. Kluberg- Stern, A. Morel, O. Napoly and B. Pctcrsson, Nucl. Phys. B 220, 447 (1983). 

D. Daniel and T.D. Kieu, Phys. Lett. B 175, 73 (1986). 



- 20 - 



